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SYSTEM AND METHOD FOR ENABLING SIMULTANEOUS CALIBRATION 
AND IMAGING OF A MEDIUM 

CROSS-REFERENCE TO A RELATED APPLICATION 
[0001] This application claims priority from the Provisional Appln. No. 

60/261,492 filed on January 12, 2001, the entire disclosure of which is incorporated 
herein by reference. 

STATEMENT AS TO FEDERALLY SPONSORED RESEARCH 
[0002] This invention was made with Government support under grants from the 

National Institutes of Health, NIH R29-NS38842 and NIH P41-RR14075, and the U.S. 
Army, No. DAMD 17-99-2-9001. As a result, the Government of the United States of 
America has certain rights to the invention, 

TECHNICAL FIELD 
[0003] The invention relates to imaging, and more particularly to calibration of 

optical methods of imaging highly scattering media. 

BACKGROUND 

[0004] Opaque media such as paint, milk, foam, emulsions, colloidal suspensions 

and human tissue do not strongly absorb visible light, while being optically turbid. This 
turbidity is a result of a very short scattering length characteristic of photons traveling 
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within these media. The light does nOt travd in straight lines through such substances, 
and instead "diffuses" in a manner similar to heat flow. In other words, photons are 
multiply scattered within these media until they are either absorbed, or they exit the 
boundaries of the medium. 

[0005] Recently, there has been a significant interest in the use of optical 

radiation for imaging within the highly scattering media, such as a biological tissue. 
Photons can travel within a highly scattering medium along a distribution of paths, of 
which very few are straight. Thus, the direction of the light into a highly scattering 
medium and a subsequent durection of the diffuse light emitted from the medium provides 
certam information regarding local variations in the scattering and absorption 
coefficients. Such information can identify, for example, a breast or brain tumor, 
bleeding in the brain, or aneurysm. Furthermore, multiple wavelsngthS Can bC USCd tO 
spectroscopically determine local tissue concentrations of oxy-hemoglobin (HbO) and 
deoxy-hemoglobin (Hb) in tissue, which can vary in response to some stunulus, e.g., a 
drug. For a general description of such apphcations, see, e.g., A. Yodh et al.. Physics 
Today, 34-40, March 1995. 

[0006] If the spatially varying optical properties of a highly scattering medium 

are known, photon propagation within the medium can be calculated numerically. This 
numerical calculation is simpUfied when the scattering is much larger than the 
absorption, in which case photon propagation can be approximated as a diffusion process. 
This condition is typically satisfied in a biological tissue in the spectral range of about 
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700 to 900 ran. The numerical calculation gives the distribution of light inside a highly 
scattering medium, and is usually referred to as the "forward calculation." For a medium 
being imaged, however, the "inverse calculation" should be solvsd, e.g., by dedUCillg ttlG 
distribution of optical properties within the medium from the difflise light measurements. 

Numerical techniques for performing the inverse Calculation Include Singular value 

decomposition (SVD), simultaneous iterative reconstruction technique (SIRT), K-space 

diffraction tomography, and a use of an extended Kaman filter. For a general review of 

techniques for the forward and inverse calculations, see, e.g., S.R. Arridge, Inverse 
Problems, \5-Ml, 1999. 

[0007] hi the art of diffuse optical tomography ("DOT"), multiple sources 

sequentially direct the light into a highly Scattering mcdium (e.g., tissue), at spatially 

separated locations. For each such source, multiple detectors on the tissue measure the 
diffuse hght emitted from the sample at spatially separated locations. The detectors may 
further obtain one or more parameters of the diffuse light emitted, e.g., fluence, and then 
utilize those parameters as input in the inverse calculation. However, the measurements 
can include various errors caused by, for example, source and detector coupling to the 
tissue, source and detector positional uncertainties, fluctuations in the source power, and 
variations in the detector gain. 

[0008] to minimize these uncertainties, DOT systems are typically calibrated with 

initial measurements for a known sample, and the calibration is used to correct the 
subsequent measurements for imaging an unknown sample. Unfortunately, errors can 
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vary from a measurement to another measurement because of, e.g., the movement or 
perspiration of a patient or the movement of an optical fiber that forms part of a source or 

detector. Thus, the results of an inverse calculation cafi include experimental Systematic 

errors caused by measurement variations that are independent of the medium's properties 
of interest. The experimental systematic errors can also limit absolute spectroscopic 
measurements of the optical properties at a particular spatial location, i.e., absolute, rather 
than relative, values of absorption and scattering. 

[0009] International Apphcation No, WO 01/1 924 1 describes a calibration 

methodology for the difiuse optical measurements that corrects the transmittance 
measurements between a source and a detector for factors unrelated to sample properties. 
The calibration methodology is based on the same set of transmittance measurements that 
are subsequently corrected by the calibration, and then are used in imaging and/or 
spectroscopy applications. This calibration method involves a forward calculation for 
each of multiple source-detector pairs based on an approximate model of the sample, and 
a minimization of an expression that depends on the forward calculations and the 
transmittance measurements to determine self-consistent coupling coefficients foi eVCiy 
source-detector pair. Once the coupling coefficients are determined, they are used to 
correct the transmittance measurements. An inverse calculation is performed on the 
corrected sample measurements to determine spatial variations in the optical properties of 
the sample. 
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SUMMARY OF THE INVENTION 
[0010] The present invention is preferably based on the fact that by allowing 

source and detector calibration factors to vary jfreely in the inverse calculation and be 
determined simultaneously with image reconstruction, image artifacts can be significantly 

reduced. 

[0011] In general, the present invention provides systems and methods for 

determining a distribution of one or more optical properties of a medium illuminated with 
radiation (e.g., electromagnetic radiation, infrared radiation, continuous-waive radiation 
etc.) from one or more sources, where the radiation exiting the medium is received by 
one or more detectors. The detectors may be adapted to obtain one or more parameters of 
the received radiation. Once the radiation is received and the parameters are obtained, 

one or more optical properties of the medium can be derived using at least one calibration 
factor and the received radiation, in which the calibration factors are variables. For a 
highly-scattering medium, the probabihty that photons entering the medium will scatter 

can significantly exceed the probability that photons entering the medium will be 
absorbed. Hence, employing diffusion approximation may be appropriate because the 
probability of scattering is at least an order of magnitude higher than the probability of 
absorption. Each source can be spatially separated firom a respective detector, e.g. by two 
centimeters, and the separation selected can be limited by the experimental parameters 
(e.g., the number of sources and detectors, the radiation used, the size of the medium 
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etc.). Where the medium is human tissue, the separation may be one or more 
centimeters. 

[0012] The calibration factors used as variables in deriving one or more optical 

properties of the medium can include, for example, a source coupling factOf, (wMch Cail 
be referred to interchangeably with a source strength or a detector gain), a detector 
coupling factor, a detector gain, a source location factor, and a detector location factor. 
[0013] The optical properties may be derived as a continuous-wave data and 

include, e.g., the absorption md scattering. However, the optical properties may also be 
derived, for example, as a frequency-domain or time-domain data. In case that the 
derived optical properties are in the frequency domain, each would have a respective 
amplitude and phase. Alternatively, if the derived optical properties are in the time 
domain, each would have a respective amplitude and temporal off-set. 
[0014] The or more of these optical properties can be derived by solving an 

inverse problem. The inverse problem can be solved by using, for example, a linear 
approximation model. One such exemplary linear approximation model used to derive 
the optical properties of the medium is a Bom Approximation. Another exemplary linear 
approximation model is a Rytov Approximation. In the Rytov Approximation model, the 
measured parameter (e.g., fluence), is defined as an exponential of a perturbation on the 
background fluence. A specific implementation of an embodiment of the method of the 
present invention using the Rytov approximation involves minimizing a least-squares 
expression for the difference between the logarithms of the theoretical and measured 
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fluence for each source-detector pair. The model can also scale the obtained parameters 
taken as input (e.g. fluence), to make those parameters dimensionless and to be on the 
same order as the calibration factors. 

[0015] Some examples of the property distributions obtained from performing the 

calculations on these measured parameters may include an absorption distribution, a 
scattering distribution, a distribution of both, etc. The determined property distribution 
can be displayed in one or more images. 

[0016] In another embodiment of the present invention, a computer-readable 

medium containing a program is provided that can cause a processor to receive radiation 
scattered through an illuminated medium as input. Also, the processor may be 
programmed so as to obtain one or more parameters of the received radiation as input. 
The program can further cause the processor to derive one or more optical properties of 
the medium by using the input and one or more calibration factors. The program can 

subsequently cause the processor to determine a property distribution using the derived 

optical properties. The property distribution may be provided as output. The program 
can use a non-linear or linear model to derive the optical properties of the media. Using 
the Rytov approximation model, the least-squares expression for the difference between 

the loganthms of thfi theoretical and measured fluence for each source-detector pair can 

be minimized so as to derive the optical properties, and determine the property 
distribution. The program can also scale all obtained parameters taken as input to make 
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tiiein dimensionless and to be of the same order as the freely varying source and detector 
calibration factors. 

[0017] In yet another embodiment of the present invention, systems are provided 

for determining the distribution of one or more optical properties of a medium 
illrnninated with radiation from one or more sources. These systems have one or more 
radiation detectors for receiving radiation. These systems also include a processor for 
deriving one or more optical properties from the received radiation and one or more 
calibration factors. The processor also determines the distribution using one or more 
properties of tiie medium. 

[0018] The sources can have optical fibers and lasers, and the detectors may 

include optical fibers linked to photodetectors. To irradiate the medium, the Systeni can 
use electromagnetic radiation, e.g., infrared radiation, or any continuous-wave radiation. 
The sources providing the radiation can te spatially separated frOffl the dCtCCtOrS that 
receive the radiation. The spatially separated sources can extend in a plane opposite to a 
plane containing the detectors, and the spatial separation can be, e.g., two centimeters, 
and can be varied according to the experimental parameters (i.e., the number of sources 
and detectors, the radifltion USed, and the size of the medium). When the medium is 
human tissue, the separation can be one or more centimeters. The processor can also 
provide absorption coefficients, scattering coefficients, or both. 
[0019] Jn addition, the property distribution determined by the system of the 

present invention can be provided as an image on a display, e.g., a computer screen. The 
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model used to determine the property distribution Can again be non-linear or linear. One 
such model includes the Rytov approximation and can involve, e.g., minimizing the 
square of the difference between the logarithms of the theoretical and measured fluence 
for each source-detector pair. The model can also utilize scaling parameters obtained 
from the signal as input to make them dimensionless and to be of the same order as the 
calibration factors. 

[0020] in still another embodiment of the present invention, systems are provided 

for determining a distribution of one or more optical properties of a medium illuminated 
with radiation from one or more sources. These systems have means for receiving 
radiation exiting the medium, and means for deriving one or more optical pfoperties Of 

the medium using one or more calibration factors and the received radiation, in which the 
calibration factors are variables. Means for determining a distribution in the medium 
using the derived optical properties are also provided. 

[0021] Li yet another embodiment of the present invention, software systems are 

provided which, when executed on a processing device, may determine a distribution of 
one or more optical properties in a medium illuminated with radiation from one or more 
sources. These systems have a processing subsystem which, when executed on the 
processing device, configures the processing device to obtain parameters of received 
radiation exiting the medium, derive one or more optical properties of the medium using 
one or more calibration factors and obtained parameters, in which the calibration factors 
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are variables, and determine a distribution in the medium using one or more derived 
optical properties. 

[0022] Unless otherwise defined, all technical and scientific terms used herein 
have the same, or substantially similar, meaning as commonly understood by one of 

ordinary skill in the art to which the present invention belongs. Although methods and 
materials similar or equivalent to those described herein can be used in the practice or 

testing of the present Invention, exemplary methods and materials are described below. 

In addition, the materials, methods, and examples are illustrative only and are in no way 

limiting. All cited references are incorporated herein by reference. 

[0023] In prior art systems and methods, several potential experimental 

systematic errors can hamper the ability to obtain images of the distribution of properties 
within a medium by applying radiation. While the characteristics of a system can be 
modeled prior to experimentation, slight fluctuations during experimentation, including 
for example changes in source strength or detector gain, can cause significant image 
artifacts. By including calibration in the image reconstruction algorithm, the present 
invention provides a substantial image improvement, despite the presence of possible 
considerable experimental systematic errors. One advantage of the methods and systems 
of the present invention is their ability to obtain diffuse optical tomography images 
having considerably reduced artifacts. The new imaging methods of the present 

invention can also be apphed in electrical impedance tomography. These imaging 
methods can be utilized to image spatially varying optical properties of biological media, 
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e.g., tissue, with enhanced contrast. In addition, such novel methods can be appHed to 
image non-biological media (e.g., plastics, ceramics, or liquids) to detect defects or 
impurities using radiation. Other features and advantages of the invention will be 
apparent from the following detailed description, and from the claims. According to the 
present invention, the source and detector calibration factors can be allowed to vary 
freely in the inverse calculation and be determined simultaneously with image 
reconstruction, and thus image artifacts can be significantly reduced. 

BRIEF DESCRIPTION OF THE DRAWINGS 
[0024] FIG. 1 is a top-level block diagram illustrating an exemplary embodiment 

of a method according to the present invention. 

[0025] FIG. 2 is a more detailed block diagram illusfrating a step of deriving 

optical properties of the method shown in FIG. 1. 

[0026] FIG. 3 is a schematic diagram of an exemplary embodiment of a diffuse 

optical tomography system according to the present invention. 

[0027] FIG. 4 is a schematic diagram of an exemplary source-detector pair for the 

system of FIG 1. 

[0028] FIG. 5 is a diagram of a first example of an experimental geometry 

derived using the present invention. 

[0029] FIGS. 6 A, 6B and 6C are exemplary images of simulated data obtained 

from the experimental geometry in FIG. 5 generated for 0%, 40% and 80% uncertainties, 
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respectively, using a normalized Rytov approximation and excluding the reconstruction 
of freely varying source and detector coupling factors. 

[0030] FIGS. 6D, 6E and 6F are exemplary images of simulated data obtained 

from the experimental geometry in FIG. 5 generated for 0%, 40% and 80% uncertainties, 
respectively, using a normalized Rytov approximation and including the reconstruction of 
freely varying source and detector coupling factors. 

[0031] FIGS. 7A, 7B, 7C and 7D are images of simulated data obtained from the 

experimental geometry in FIG. 5 generated for 0%, 40%, 80% and 120% uncertaiuties, 
respectively, using a Bom approximation and including the reconstruction of freely 
varying source and detector coupling factors. 

[0032] FIGS. 7E and 7F are images of simulated data obtained from the 

experimental geometry in FIG. 5 generated for 80% and 120% uncertainties, respectively, 
using a normalized Rytov approximation and including the reconstruction of freely 
varying source and detector coupling factors. 

[0033] FIGS. 8A, 8B, 8C, and 8D are images, similar to those of FIGS 7A-7F, of 

simulated data obtained from the experimental geometry in FIG. 5 generated for 10%, 

20%, 40% and 60% uncertainties, respectively, using a R540V approximation with the 
reconstruction of freely varying source and detector coupling factors, but without the 
normalization. 

[0034] FIG. 9A is a set of images of simulated data obtained from the 

experimental geometry in FIG. 5, generated for localization uncertainty of e.g., 0.1 mm 
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using a normalized Rytov approximation and excluding the reconstruction of freely 
varying source and detector location factors. 

[0035] FIG. 9B is a set of images of simulated data obtained from the 

experimental geometry in FIG. 5, generated for localization uncertainty Of e.g., 0.1 IHIIl 
using a normalized Rytov approximation and including the reconstruction of freely 
varying source and detector location factors. 

[0036] FIG. 1 OA is a set of images of simulated data obtained from the 

experimental geometry in FIG. 5, generated for localization uncertainty of e.g., 0.5 mm 
using a normalized Rytov approximation and excluding the reconstruction of freely 

varying source and detector location factors. 

[0037] FIG. 1 OB is a set of images of simulated data obtained from the 

experimental geometry in FIG. 5, generated for localization uncertainty of e.g., 0.5 mm 
using a normalized Rytov approximation and including the reconstruction of freely 
varying source and detector location factors. 

[0038] FIG. 1 1 A is a set of images of simulated data obtained from the 

experimental geometry in FIG. 5, generated for localization imcertainty of e.g., 1.0 IMl 
using a normalized Rytov approximation and excludmg the reconstruction of freely 
varying source and detectof lOCation faCtOIS. 

[0039] FIG. IIB is a set of images of simulated data obtained from the 

experimental geometry in FIG. 5, generated for locahzation uncertainty of e.g., 1.0 mm 
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using a normalized Rytov approximation and including the reconstruction of freely 
varying source and detector location factors. 

[0040] FIG. 12 is an example of the experimental setup of tke system according 

to tiie present invention. 

[0041] FIG. 1 3 A is a set of images of data obtained from the experiment of FIGS. 

9A and 9B using a normalized Rytov approximation excluding the reconstruction of 

freely vafying source and detector calibration factors. 

[0042] FIG. 13B is a set of images of data obtained from the experiment of FIGS. 

9A and 9B using a normalized Rytov approximation including the reconstruction of 
freely varying source and detector calibration factors. 

DETAILED DESCRIPTION 
[0043] The novel methods of the present invention calibrate sources and detectors 

as part of the inverse calculation for image reconstruction to reduce image artifacts. 

General Methodology 

[0044] Research in diffiise optical tomography has rapidly progressed in moving 

the imaging modahty from the development stage of computer simulations and phantom 
experiments to the application stage of imaging animal and human subjects. However, 
little work has discussed the importance and solution of experimental details that degrade 
image quality including, without limitation, the freatment of boundary conditions for 
solving the forward problem, source and detector coupling to tissue, source and detector 
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positional uncertainties, and intrinsic tissue heterogeneity. The novel calibration 
techniques according to the present invention are based on the principle that, e.g., the 
inclusion of calibration factors as variables in the image reconstruction algorithm 
significantly improves the image quality. 

[0045] The calibration techniques of the present invention can be applied to 

determine a property distribution within a scattering medium, which is an ill-posed, non- 
linear inverse problem. A linear approximation can be obtained by assuming the 
v^ations may be described as small perturbations in absorption and scattering, Spa and 
5]is\ respectively, around known background values of absorption and scattering, Uao 
and Uso', respectively. As indicated previously, one or more parameters of diffiise light 
emitted through a scattering medium may be obtained and used as input in the inverse 
calculation. One of those parameters is preferably fluence. A linear approximation of the 

ffleasured fluence 0 can be obtained as a perturbation on a background fluence %. 

In imaging, the fluence 0 is measured experimentally and the perturbation 4>pert (which is 
caused by spatial variation in the radiative properties of the sample) should preferably be 
calculated. 

[0046] An accurate determination of the perturbation ^p^n prefers a calibration of 

a source and a detector, and a reasonable estimate of the background fluence ^o- Errors 
in the determination of the calibration factors are likely to result in a highly localized 
absorption and scattering perturbations appearing adjacent to the sources and detectors as 
a compensation for the errors in the image reconstruction process. This type of image 
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noise can be evidenced by high frequency spikes appearing preferentially near source and 
detector locations. 

[0047] The exemplary methods of the present invention allow for the calibration 

by enablmg the calibration factors to vary freely in the inverse calculation, and 
determining them as part of the image reconstruction algorithm. There are various 
calibration factors, e.g., source and detector coupling factors j and d, respectively, and 
source and detector location factors rj and fd , respectively. By considering the logarithm 
of the measured fluence <l> as input into the image reconstruction algorithm, the 
perturbation Open becomes linearly dependent on the logarithms of s and d. However, the 
perturbation Opert can remain non-linearly dependent on the optical properties of the 
sample. To obtain a linear dependence on the optical properties, an approximation, e.g., 
the Rytov or Bom approximation, can be used, and the problem thus can become 
completely linear. 

The Forward Problem: 

[0048] The radiative transport equation provides a rigorous theory to describe 

radiation emitted from a medium. It may be applied to the migration of photons through 
highly scattering media, such as biological tissue. This approach can indicate that near- 
infrared photons in highly scattering media e.g., human and animal tissue, essentially 
undergo a random walk, since the scattering probability greatly exceeds the absorption 
probability. The propagation of these photons through a highly scattering medium can 
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therefore be described by a diffusion approximation to the radiative transport equation. 
This exemplary diffusion can be provided as an equation as follows: 

-V-(D(r)VO(r,t)) + v)a rr)0(r,t)+-^^^- vS(r,t), (1) 
dt 

where D(r) is the photon diffusion coefficient, 6(r,t) is the photOn flueilCe at pOSltiOIl I 
and time t, v is the speed of light in the medium, Ua(r) is the absorption coefficient, and 
S(Y,t) is the source diStriljUtion of photons (see, e.g., A. Ishimaru, Wave Propagation and 
Scattering in Random Media, Academic Press, Inc. San Diego, 1978; M. S. Patterson et 

al., Applied Optics 28:2331, 1989; R. C. Haskell et al, Journal of the Optical Society of 
America A, \\\T1H, 1994). The photon diffusion coefficient D{r) is defined by the 
equation: 

D{r)= / — , (2) 

3^i/(r)+a 

where ]Xsi^) is the reduced scattering coefficient. The coefficient a can be variously set 
to three, zero, or some other value (see, e.g., D. J. Durian, Optics Letters, 25:1502, 1Q98; 
R. Aronson and N, Comgold, Journal of the Optical Society of America A, 16:1066, 
1999). 

[0049] When focusing on the variations in the absorption, the value selected for a 

Is moot, and therefore it can preferably be set to zero. 

[0050] When the distribution of properties within the medium varies spatially, 

two approaches can be used to find approximate solutions to Equation (1): the Bom 
approximation (see, e.g., A. Ishimaru, Wave Propagation and Scattering in Random 
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Media, Academic Press, Inc. San Diego, 1978; A. C. Kak and M. Slaney, Principles of 
Computerized Tomographic Imaging, IEEE Press, New York, 1988), and the Rytov 
approximation (see, e.g., A. C. Kak and M. Slaney, Principles of Computerized 
Tomographic Imaging, IEEE Press, New York, 1988). The Bom approximation can 
preferably be defined as: 

The Rytov approximation is preferably : 

<1)=^„ exp((I)p^rt)- (4) 
[0051] Each approach decomposes the total fluence ^ into the background 

fluence <t>o, which only depends on the background optical properties and -\Xso\ and 
(l>pert, which is linearly related to the spatial variations in the optical properties Sjia and 
5Us. For the Rytov approximation, assuming absorption variations only, the perturbation 
can be calculated as follows: 

|<I..(r„r)i«|«(;(r,r,)*. (5) 

where and are the respective positions of the source and detector, Do is the diffusion 
coefficient for the chosen a, and G is the Green's function of the photon diffusion 
equation for the background optical properties provided by the boundary conditions. If 
the background is homogeneous, G can be expressed analytically for some simple 
geometries (see, e.g., M. S. Patterson et al.. Applied Optics 28:2331, 1989; R C. Haskell 
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et al., Journal of the Optical Society of America A, \ l:21TI, 1994; S. R. Arridge et al., 
Physics in Medicine and Biology, 37:1531, 1992). 

[0052] One exemplary geometry to which Equation (4) can be applied is a slab, 

i.e., a geometry where the sources are in a single plane on one side of the sample and the 
detectors are in a plane on the opposite side thereof. In this geometry, the extrapolated 
zero-boundary condition can be applied to the above-described radiative transport 
problem to solve for the Green's fimction G of the diffusion equation for the background 
optical properties (see, e.g., R. C. Haskell et al., Journal of the Optical Society of 
America A, 11:2727, 1994). The background fluence i>o IS related tO the Giecn'S ftUlCtlon 
by the source and detector coupling factors. The source coupling factor s includes source 
power and a coupling to the medium. The detector coupling factor d includes detector 
gain and a coupHng to the medium. Using 5 and d respectively, to represent the source 
and detector coupling factors, the background fluence is as follows: 

<^,(r„r,)=sdG{r„r,) (6) 

The Inverse Problem: 

[0053] By solving the inverse problem, the images of spatially varying optical 

properties can be generated from measurements of the fluence The inverse problem 
includes solving a least-squares equation. For the Rytov approximation, the expression to 
minimize can be as follows: 

F(x) = 5[ln O^^^, (x) - In^,^. } , (7) 
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where the index / is summed over the measurements for each source-detector pair, (^i(x) 

is provided by Equations (4) and (5) where x is a vector giving 5]ia at each voxel 

position, and <J>Meas(x) is the measured fluence for each source-detector pair. 

[0054] For the case of fewer measurements than unknowns, the linear problem is 

underdetermined, and can be described by the (regularized) Moore-Penrose generalized 

inverse: 

x=-A^{aA^ +llY y, (8) 
where / is the identity matrix, and X is the Tikhonov regularization parameter. Each 
element of the matrix A can be defined as: 

where and rd,i are the positions of the i* source and the i* detector, respectively. The 
position of the j^^ VOXel is eXpreSSCd as Yj. Using the Rytov approximation, each element 
of ;v can be as follows: 



J; =ln 



(10) 



[0055] In one exemplary approach, X can be set to 1 0"^ of the maximum eigen- 

value of^^(see, e.g., H. Dehghani et al.. Physiological Measurement, 20:87, 1999; Y. 
Kolehmainen et al., Physiological Measurement, 18:289, 1997). 
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Image and Calibration Factor Reconstruction 

CALIBRATION FOR THE SOURCE AND DETECTOR 
COUPLING FACTORS AND MAGE RECONSTRUCTION 

[00561 Image quality depends on accurate knowledge of the backgroimd optical 

properties, }iao and }iso', and the source and detector calibration factors s and d. 
Therefore, estimates of the matrices A and y similarly depend Oil thCSC quantities. By 
considering and ]iso' as known properties, the inverse problem can be solved for the 

calibration factors. 

[00571 As described above, the measurements can include various errors caused 

by, for example, source and detector coupling to the medium, fluctuations in source 
power, variations in detector gain, and source and detector positional uncertainties. 
Hence, the measurements can be calibrated for the source and detector coupling and 
intensity fluctuations (source and detector coupling factors), as well as for their positions 
( source and detector location factors). 

[0058] In the case when measurements should be calibrated for the source and 

detector coupling factors s and d (the source and detector coupHng and fluctuations in tiie 
source power and the detector gain), an augmented model of a particular measurement, 
may be written in terms of the vuiknown source coupling factors Sk, the detector coupling 
factors di, and the absorption perturbations Sucy. The source and detector coupUng 
factors can preferably refer to both coupling fluctuations and intensity fluctuations of the 
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sources and detectors. Using these variables, the particular measured parameter (in this 

case fluence) may be expressed as follows; 



=ln 



(11) 



where k(i) and identify the source and detector, respectively, which are used for the 

measurement, 

[0059] In case that numerous measurements are conducted, the measured 

parameter (fluence) can be written in a matrix form as: 

y=B^,wh6r6B=[ASD], (12) 

where S is a source coupling factor matrix, D is a detector coupling factor matrix, and A 
is a matrix A rescaled by a factor of jiao. 
Factor ^ is defined by the formula: 



• Insi ••• Inj^^ •-• hn^i ■■■ ln</^^ 



(13) 



where A^v is the number of voxels, Ns is the number of sources, and Nd is the number of 
detectors. Sucy- can be scaled by }iao to make the elements dimensionless and being of the 
same order as In 5 and hi d. ScaUng these quantities to make them dimensionless and of 
the same order as the logarithms of the calibration parameters yields a ttOftHalized 
formulation. 

[0060] As indicated previously, A is a rescaling of the matrix A whose terms are 

defined by Equation (9) by a factor of Uao. Thus, 



NY02:364961.1 



-22- 



A34927- 069225.0110 



A=^t„„A. (14) 

[0061] Matrices S and D have block diagonal form. For example, S haS a OflB ill 

the column for each measurement corresponding to source k, and D has a one in the 
column for each measurement corresponding tO deteCtOI /. ThUS, With fOUT measurements 
between two sources and two detectors, the following matrix is defined: 



[SD]= 



10 10 

10 0 1 
Olio 
0 10 1 



CALIBRATION FOR THE SOURCE AND DETECTOR LOCATION 
FACTORS AND MAGE RECONSTRUCTION 

[0062] Conversely, if the measurements should be calibrated for the source and 

detector locations, an augmented model of a particular measurement, >>„ may be defined 
in terms of the uncertain source location 5rsk, uncertain detector location 5rdk and the 
absorption perturbations Suaj- In this exemplary augmented model, yt may be expressed 



yi =in 



+ x-r ^V,,(,)a)(r,,,(.,,r,,(J.8r,,(,) + ^ Aj^i'aj^ 



(15) 



where k(i) and l(t) identify the source and detector, respectively, used for the / 
measurement. It should be noted that 5rsk and 6rdi represent the differences from the 
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respective assumed values as used in and the real values for the measurement of <l>. In 
addition, 5rsk and br^ can be vector quantities representing the uncertainties in the x, y, 
and Z coordinates of the corresponding source or detector. 

[0063] For this exemplary embodiment of the present invention, in case that 

numerous measurements are conducted, the measured parameter (fluence) can be written 
in a matrix form as follows: 



y=B^, whereB=[APsPd] 
where E, is defined by: 



5r, 



(16) 



(17) 



in which Nv is the number of voxels, is the number of sources, and Nd is the number of 
detectors. Scaling b]iaj by ]iao and the positional unknowns by cr makes the elements 
dimensionless. a can be selected so as to make the unknowns of the same order of 
magnitude, and in this manner, a is preferably an estimate of the expected error in the 
source and detector locations. Similarly, A is a rescaling of the matrix A whose terms are 
given by Equation (9) by a factor of vioto- Thus, 

A=Ma„A. (18) 

[00641 Matrices Ps and Pd can have a block diagonal form. For example, Pg may 

be non-zero in the column for each measurement corresponding to source k, and Pd 
may be non-zero in the f ^ column for each measurement corresponding to detector /. 
These preferably have the same form as the [S D] matrix provided above with respect to 
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the image reconstruction including the caUbration for the source and detector coupling 
factors. 

Frequency Domain Formulation 

[0065] The imaging method according to the present invention is not limited to 

using continuous-wave radiation, hi the frequency domain, the measured parameters of 
amplitude and phase may also be obtained. If the measured parameters are in the 
frequency domain, then, for example, the source and detector coupling factors have 
corresponding unknown phases 6^ and 6^, respectively. They also likely have unknown 
amplitudes As and Ad, The in s and In d, which appear in Equation (13), can be complex, 
as are the corresponding elements in the matrix A. To separate the real and imaginary 
portions of the complex matrix, the real portions can be stacked on top of the unaginary 
portions, i.e., 



re (a) re{s) re[D) 

imU) imiS) im{D) 



(19) 



[0066] hi the frequency domain, should be adjusted as follows: 



1 = 



(20) 
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[0067] For the case of four measurements between two sources and two detectors, 

the complex matrix can be: 



[SD]= 



10/010/0 

1 0 / 0 1 0 I 0 

0 1 0 z 0 1 0 z 

0 1 0 i 0 1 0 I 



m) 



[00681 The example of a formulation for the case wheffi there are Uncertain 

source and detector locations si substantially similar to the example provided above. In 
particular, Equation (19) would be provided as follows: 
refe) re(Ps) re(pj' 



im(A) im(Ps) im(Pi,) 



(22) 



Time Domain Formulation 

[0069] The imaging method according to the present invention can also be 

applied to measurements taken in the time domain. To obtain time domain data, each 
source provides a temporally coherent pulse of radiation, e.g., a picosecond light pulse, 

and the detectors are time-gated to measure the temporal delay of the radiation in addition 

to its intensity. For a general reference on such time-domain techniques in diffuse optical 
tomography see, e.g., M.S. PatterfiOll et al.. Applied OpUCS 28:2331, 19§9; S.R. Arridge, 
Inverse Problems, 15:841 , 1999. Therefore, to apply the new calibration methods to 
time-domain data, the calibration factors s and d are described by corresponding 
unknown amplitudes and temporal off-sets. The unknown amplitudes in the time-domain 
are treated in the same fashion as in the continuous- wave domain. The unknown 
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temporal off-sets represent a non-linear perturbation and can be reconstructed using a 
non-linear form of the model. 

[0070] The example of a formulation for the case where there are uncertain 

source ^d detector locations si substantially similar to the example provided above. 

Exemplary Illustrations of the System and Method of the Present Invention 
[0071] Referring to FIG. 1 , a top-level block diagram presenting an exemplary 

embodiment of a method according to the present invention is illustrated. Initially, the 
radiation which exits from the medium is received (step 110). Then, in step 120, one or 
more optical properties of the medium are derived using the radiation received in step 
110 and one or more calibration factors, wherein the calibration factors are variables. 
Subsequently, in step 130, a distribution of the optical properties in the medium is 
determined using the optical properties dfilived In Step 120. ThC details Of Steps 110- 
130 are described in detail herein above. 

[0072] Further details for step 120 of deriving one or more optical properties of 

the medium are additionally illustrated in FIG. 2 and described below. In particular, 
parameters of radiation are obtained in step 210. One of such parameters can be fluence. 
The detectors with which the radiation exiting the medium is received maybe adapted to 
directly measure the parameters of the radiation, e.g., fluence. Alternatively, the 
parameters maybe calculated using the data corresponding to the received radiation. 
Once the parameters are obtained in step 210, the formulation of each obtained parameter 
can be augmented to include the calibration factors as variables (step 220). As described 
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above, the calibration factors may include the source coupling factor, detector coupling 
factor, source location factor and detector location faCtOT. SubSCqUSIlt tO augnenting the 
formulations of the obtained parameters of radiation, one or more optical properties of the 
medium are derived as output by SOlving an inveise problem in which the augmented 
formulations are used as input (step 230). 

General System 

[0073] FIG. 3 shows an exemplary embodiment of the system capable of 

executing the method according to the present invention in which experimental 
systematic errors may cause significSflt image artifacts. SUCh System Can be a diffuse 
optical tomography ("DOT") system 300. The system 300 includes an array of spatially 
separated light SOUfCes 310 and Spatially separated detectors 320. During use, the array 
of somrces 310 and detectors 320 can be positioned over a sample 350 to be imaged, e.g., 
a patient's head or breast. A controller 330 connected to the light sources 310 can 
sequentially trigger them to forward or provide light into the sample 350, so that a highly 
scattering medium may cause the light to become di^^se within the sattipls. FOI each 
sequentially triggered source, each detector 320 receives hght that reaches it through 
sample 350 to obtain the measured fluence The controller 330 is also connected to the 
detectors 320, and selectively channels the signals from the detectors 320. An analyzer 
340 is connected to the controller 330 and analyzes the signals generated by detectors 
320. 

-28- 
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[0074] For example, in one exemplary embodiment shown in FIG. 4, each source 

310 may include a diode laser 400 for producing optical radiation and a lens 420 for 
coupling the optical radiation into an optical fiber 410. The optical fiber 410 has an end 
415 adjacent to the sample 450 for directing the optical radiation into the sample 450. 
Each detector 420 includes an Optical fiber 452 having an end 455 adjacent to the sample 
450 for receiving the optical radiation exiting the sample 450, a photodetector 460 for 
measuring the intensity of the optical radiation received by the optical fiber 452, and an 
amplifier 470 for amplifying the output of the photodetector 460. 
[0075] To obtain information about the optical properties of the sample 450 firom 

the fluence, it is preferable to consider the measured fluence 0 as representing the 

perturbation ^pen on a background fluence Oo. According to Equation (4), this can be 
accomphshed by using the Rytov approximation. By minimizing the least-squares 
expression as provided in Equation (7), the analyzer 340 can simultaneously obtain the 
optical properties of the sample and conduct source and detector caUbration. 
[0076] In other embodiments of the present invention, the light source can include 

a laser other than a diode laser, e.g., an ultrafast laser, or may be an incoherent source. 
Also, the sources can include a common light source that selectively couples light into 
multiple fibers that deliver the light to spatially separated locations on the sample. 
Alternatively, the sources need not include optical fibers at aU. For example, the lasers 
themselves can be positioned adjacent to the sample or can include beam delivery optics 
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to direct the light to the sample tlirough free space. Furthermore, the hght sources can 
provide light at multiple wavelengths by including, for example, multiple diode lasers. 
[0077] In the description provided below, it can be assumed that the sources 

provide continuous-wave optical radiation, and that the detectors measure tiie intensity of 
the optical radiation. 

[0078] It should be understood that the calibration techniques described herein 

can also be appHed to other measurement techniques in which the sources do not provide 
continuous-wave radiation. For example, in some techniques, the amplitude of optical 
radiation provided by the sources can be modulated to create photon density waves in the 
sample, and the detectors are configured to measure the amplitude and phase of the 
photon density waves after propagation through the sample. For a general reference on 
"DOT" techniques with modulated optical radiation, see, e.g., M.A. O'Leary et al., Phys. 
Rev. Lett. 69, 2658 (1992). Furthermore, in other techniques, each source can provide a 
temporally coherent light pulse, for example, a picosecond pulse, and the detectors may 
be time-gated to measure the temporal delay of the diffuse light pulse in addition to its 
intensity. For a general reference on such time-domain DOT techniques, see, e.g., M.S. 
Patterson et al, Appl. Opt. 28:2331, 1989, and S.R. Arridge, Inverse Problems, 15:841, 
1999. 

Exemplary Applications 

[0079] The methods and systems of the present invention provide the ability to 

obtain images of spatially varying optical properties within highly scattering media, such 
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as human and ailifflal tisSUe, and can therefore be used to address many new biomedical 
problems and applications. For example, these novel methods and systems can be used to 
image a region of tissue to ascertain the presence of an early tumor, a small amount of 
bleeding, or an early aneurysm. In another exemplary application, multiple wavelengths 
can be used to evaluate local tissue concentrations of chemicals, such as hemoglobin, 
within tissue. Also, the methods and systems of the present invention can be fiirther 
combined with use of a chemical stimulus to determine the chemical response, 
[0080] Non-biomedical applications for the methods and systems of the present 
invention also exist and can be used to image materials such as plastics, ceramics, and 
hquids for defects or impurities. In addition, the methods and systems of the present 
invention can be used for identifying regions of interest having distinct optical properties 
within highly scattering media, such as concealed objects, without invasive study of that 
media. For example, the methods of the present invention can be used tO image an 
optical heterogeneity within a diffusive liquid. Furthermore, the novel calibration 
methods and systems can also be applied in electrical impedance tomography. 

Software Implementation 

[0081] The calibration methods can be implemented in hardware, software, or a 

combination of both. The methods can be implemented via computer programs using 
standard programming techniques following the methods, equations, and Figures 
described herein. Program code can be executed using a processing arrangement which 
can receive input data to perform the ftinctions described herein and generate output 
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information. The output information can be applied to one or more output devices, such 
as a display monitor, printer etc. 

[OOS^l Each program can preferably be implemented in a high-level procedural or 

object-oriented programming language to communicate with a general computer system. 
However, the programs can be implemented in assembly or machine language, if desired. 
In any case, the computer language can be a compiled or interpreted language. 
[0083] Each computer program can preferably be stored on a storage medium or 

device (e.g., ROM, magnetic diskette, or optical disc) readable by a general or special 
purpose programmable processing arrangement, for configuring and operating the 
computer when the storage media or device is read by the processing arrangement tO 
perform the procedures and techniques described herein. The computer program can also 
reside in cache or main memory during program execution. The piOCeSSillg methodS Cao 
also be implemented as a computer-readable storage medium, configured with a 
computer program, where the Storage medium so configured can cause the processing 
arrangement to operate in a specific and predefined manner so as to perform the functions 
and techniques described herein. 

[0084] For example, referring again to FIG. 3, the analyzer 340 includes a 

processor 370, an input/output control card 360, a user interface 390 such as a keyboard 
and monitor, and memory 380. The memory 380 StoreS a program 385 specifying the 
steps of the calibration method. When executed, the program 385 causes the processor 
370 to carry out the steps of the method according to the present invention. 
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EXAMPLES 

[0085] The methods and systems of the present invention are further described m 

the following examples, which do not hmit the scope of the invention described in the 
claims and are provided merely for illustrative purposes. 

Example 1 — Simulation of Image Reconstruction 

[0086] Several advantages of the methods of the present invention were 

demonstrated by simulations. The results obtained utilizing 3 normalized RytOV 
approximation without the reconstruction of calibration factors (i.e., without using the 
new calibration methods), a normalized Rytov approximation with reconstruction, a Bom 
approximation with reconstruction, and a Rytov approximation without normalization 
were compared. A diagram of a first example of an experimental geometry is shown in 
FIG. 5. 

[0087] The computer simulation Involved the transmission through a 6-cm thick 

slab with background optical properties v^so' = 10 cm"^ and Uao, = 0.15 cin\ A 1.6-cm 
diameter absorbing object 550 with Uso' = 10 cm"^ and = 0.15 cm"^ was centered at (x, 
y, z) = (1, -1, 3) cm. Sixteen sources 510 were programmed in a four-by-four grid at z = 
0, spaced 2 cm apart to span x = -3 tO 3 CHI and y = -3 tO 3 Cm. Similarly, sixteen 

detectors 520 were located in a four-by-four grid at z = 6, spanning x = -3 to 3 cm and y 
= -3 to 3 cm in 2-cm steps. All simulated measiu^ements were made with continuous- 
wave sources. 
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[0088] The Full Bom expansion was used to solve the forward problem. A 1 .6 x 

1 .6 X 1 .6 cm cube 555 divided into 7x7x7 voxels was centered over the 1 .6-cm 
diameter absorbing object 550. 

[0089] Initially, the source and detector coupling factors were varied, namely, the 
source and detector amplitudes. The source and detector amplitudes were chosen 
randomly from a normal distribution with a mean of I . The effect of different standard 
deviations (0%, 40%, and 80%) was investigated. The normal distribution was biased by 
discarding negative amplitudes. A separate instance of the source and detector amplitude 
variation was considered for each standard deviation. There was no additive 
measurement noise (i.e., shot noise or electronic noise) in the simulated data, just the 
model error associated with the source and detector amplitudes. 
[0090] For the inverse problem, voxel centers were distributed from x = -3 to 3 

and y = -3 to 3 in 0.5-cm StepS and Z = 0.5 to 5.5 in 0.5-cm steps. There were 1859 
voxels, each with dimensions of 0.5 x 0.5 x 0.5 cm. The source and detector amplitudes 
were assimied to be all equal to 1. Two different inverse problems were considered: (1) 
Rytov without coupling, and (2) Rytov with coupling. In each case, Tikhonov 
regularization was used to obtain the pseudo-inverse, since the system matrix is under- 
determined and ill-conditioned. 

[0091] FIGS. 6A - 6F show the results obtained with the normalized Rytov 

approximation both excluding the reconstruction of the soxurce and detectors coupling 
factors (FIGS 6A-6C), and including such factors (FIGS. 6D-6F). FIGS. 6A and 6D 
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show results for 0% uncertainty in source and detector strengths, FIGS. 6B and 6E depict 
images for 40% uncertainty, and FIGS. 6C and 6F reflect 80% uncertainty. All images 

span X and Y from -3 to 3 cm, and Z-slices are vertically arranged from 0.5 (top) to 5.5 

cm (bottom). 

[0092] The results in FIGS . 6A - 6C show images without the reconstruction of 

the source and detectors coupling factors, i.e., the effect of neglecting experimental 
systematic errors in the calibration. With no variance (i.e., no systematic error in the 
sources and detectors), the object is localized to the correct voxd in X and Y and iS 
within 1 voxel (0.5 cm) in Z. The amplitude of the object is slightly underestimated due 
to blurring. While it is expected that Spa = 0.10 cm'', the reconstruction yields Siia = 
0.098 cm"^ The image contains artifacts near the sources and detectors. These artifacts 
disappear when the Bom approximation is used for the forward problem, indicating that it 
arises from mismatch between the exact forward solution and the Rytov approximation 
for the inverse model (i.e., it results from model noise unrelated to source and detector 
calibration factors). 

[0093] As the standard deviation in the calibration factors increases from FIG. 6A 

to FIG. 6C, it is still possible to ascertain the object, but this deviation has reduced the 
contrast relative to the surrounding voxels, and its amplitude is smaller than that in the 
voxels closest to the soiu-ces and detectors. The artifact amplitudes may greatly exceed 
the displayed contrast, and have been truncated to preserve contrast sensitivity sufficient 
to reveal the object of interest. The model noise caused by an incorrect calibration of the 
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source and detector calibration factors (i.e. systematic errors in s and d) clearly degrades 
the image quality by increasing the artifacts. 

[0094] FIGS. 6D - 6F show the reconstruction results where the source and 

detector coupling factors are simultaneously determined with the optical properties of the 
sample using the novel calibration methods of the present invention. Without the 
variance in the calibration factors, the location of the object remains properly resolved in 
X and Y coordinates and is within 0.5 cm in Z coordinate. The resolution in FIGS. 6D - 
6F, with the simultaneous reconstruction of the calibration factors, is poorer than in 
FIGS. 6A - 6C, due to its lacking this simultaneous reconstruction. The presence of more 
unknowns in the inverse calculation using simultaneous reconstruction explains the 
decrease In the resolution. By Contrast, the artifact shown in FIG. 6 A due to the error 
between the exact and the Rytov approximation to the forward problem is reduced using 
the simultaneous reconstruction. The reconstruction appears to compensate for the model 
mismatch in the Rytov approximation. 

[0095] Comparing the images generated using the simultaneous reconstruction of 

the source and detector coupling factors (shown in figs. 6D - 6F) to the corresponding 
images generated without reconstruction (shown in FIGS. 6 A - 6C) reveals a significantly 
improved image quality using the reconstruction. For both 40% and 80% variance, the 
image quality using the reconstruction in fact remains as good as the case where there is 
no variance. The minor differences are masked within the contrast sensitivity. While the 
reconstruction may not accurately determine the independent source and detector 
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coupling factors with the sample, it assists with the accurate determination of the source- 
detector coefficient products (e.g., Sk di) for each measurement. Despite the uncertainties 
of 80%, the results shown in FIG. 6D - 6F indicate that the source-detector coefficient 
products CSB. be determined within an accuracy of a few percent. In addition, it has been 
observed that, when additive shot or electronic measurement noise is included in the 
simulation, there may be a greater sensitivity to the measurement noise as the source and 

detector variance increases. 

[0096] FIGS. 7 A - 7F show images of simulated data obtained from the 

experimental geometry in FIG. 5 generated for 0%, 40%, 80% and 120% uncertainties, 
respectively, using a Bom approximation including the reconstruction of freely varying 
source and detector coupling factors. FIGS. 7E and 7F are images of simulated data 
obtained from the experimental geometry in FIG. 5 generated for 80% and 120% 
uncertainties, respectively, using a normalized Rytov approximation including the 
reconstruction of freely varying source and detector coupling factors. FIGS. 7A - 7F are 
usefiil to compare use of the Bom approximation with the reconstruction of the source 
and detector coupling factors to use of the normalized Rytov approximation with the 
reconstruction. FIGS. 7A - 7D show results for 0% imcertainty to 120% uncertainty in 
the source and detector strengths. By comparison, FIGS. 7E and 7F depict images for 
80% and 120% uncertainty using the Rytov approxunation. The Z-slices are vertically 
arranged with each image spanning X and Y axes from -3 to 3 cm. 
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[0097] Images generated using the Bom approximation are clearly inferior to 

those generated using the Rytov ^proximation for the same level of UHCertaillty. Meed, 
the images generated using the Rytov approximation for uncertainties of 80% (FIG. 7E) 
and 120% (FIG. 7F) resemble those generated using ilie Bom approximation at 0% (FIG. 
7A) and 40% uncertainty (FIG. 7B) more than images using the Bom approximation at a 
comparable level of uncertainty (FIGS. 7C and 7D). The standard Bom implementation 
of the image reconstmction problem is non-linear with respect to the source and detector 
coupUng factors. Therefore, if the initial estimate for the coupling factors is off by more 
than approximately 10%, simultaneously solving for the coupling factors may not lead to 
an improved result. By using the Rytov implementation (i.e., taking the logarithm of the 
data), the dependence on the soiirce and detector factors is made linear, and it is then 
possible to accurately reconstract the coupling factors despite an initial guess, possibly 
off by several hundred percent. 

[0100] While merely substituting the Rytov approximation for tke Bom 

approximation offers an improved image quality, the normalization of the different 
factors in the inverse calculation (i.e., absorption, scattering, and calibration factors) can 
also contribute to the enhanced image quality. The normaUzation has a beneficial effect 
because the reconstmction of several different factors simultaneously is difficult if the 
factors differ in magnitude by an appreciable amount. Using the normalization, the 
reconstmction can proceed accurately over a larger range of values. FIGS. 8A - 8D show 
the exemplary images using the Rytov approximation with the reconstmction of the 
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source and detector coupling factors, but without the nomialization. FIG. 8A shows the 
results for 10% uncertainty, FIG. 8B for 20%, FIG. 8C for 40%, and FIG. 8D for 60%: 
FIGS. 8A - 8D have the same X-Y scale and arrangement of Z slices as in FIGS. 7A to 
7F for an ease of comparison. Comparing FIGS. 8A - 8D with FIGS. 7E and 7F reveals 
that normalizing results in the improved image quality, even at highCI ICVClS Of 
uncertainty. 

[0101] Alternatively, a different set of calibration factors may be varied. For 

example, the source and detector locations may be varied randomly with a normal 
distribution with a mean of 0. The effect of different distribution widths (0. 1 , 0,5, and 1 .0 
mm) was investigated. A separate instance of the source and detector amphtude variation 
was considered for each distribution width. There was no additive measurement noise 
(i.e. a source or detector electronic noise) in the simulated data, and only the model error 
associated with the source and detector locations was added. 

[0102] For the inverse problem, voxel centers were distributed from x=-3 to 3 and 

y=-3 to 3 in 0.5 cm steps and z=0.5 to 5.5 in 0.5 cm steps. Approximately 1859 voxels 
were provided with dimensions of 0.5 x 0.5 x 0.5 cm. The vector ^ in Equation (14) was 
initialized to zeros. Two different inverse problems were considered: (1) the Rytov 
approximation without the localization, and (2) the Rytov approximation with the 
localization. In each case, the Tikhonov regularization was used to obtain the pseudo- 
inverse since the system matrix was imder-determined and ill-conditioned. 
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[0103] FIGS. 9A-1 IB illustrate the results obtained with the Rytov approximation 

excluding and including the reconstruction of the coupling coefficients for optode 
localization errors of 0.1, 0.5 and 1.0 mm, respectively. Referring to FIGS. 9A and 9B, 
the unages of the simulated data obtained fi-om the experimental geometry shown in FIG. 
5 usmg a normalized Rytov approximation excluding and including the reCOnStrUCtiOIl Of 
freely varying source and detector location factors, are illustrated for localization 
uncertainty of 0.1 mm. With the localization error of 0.1 mm, it is possible to locate the 
position of the absorbing heterogeneity despite the presence of image speckle near the top 

and. bottom surfaces of the scattering medium. 

[0104] With the localization errors of 0,5 mm and greater (See FIGS. 1 OA and 

lOB, and FIGS. 1 1 A and 1 IB, the surface speckle noise obscures the presence of the 
absorbing heterogeneity. The speckle noise may arise directly from the localization error. 
In most, if not all, cases, including the calibration of optode location in the inverse 
problem significantly reduces the surface speckle noise revealing the absorbing object. 

Example 2 — Phantom Results 

[0105] The actual image results of a phantom have also been obtained, and FIG. 

12 shows the experimental setup for imaging the phantom. A phantom box 1210 has a 
set of sources and a set of detectors on its top and bottom plates, respectively. The 
phantom box 1210 is connected with an RF Imager 1230 with a set of detector fibers 
1220. The phantom box 1210 is also connected with a controller 1280 via a set of source 
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fibers 1270. A processing arrangement 1290 is coupled with the controller 1280. The 
phantom box 1210 is further coupled with a beaker 1240 and a pump 1250. 
[0106] The phantom box 1210 was filled with 0.5% intralipid solution, and a 2- 

cm diameter heterogeneity sphere was embedded in this solution. The content of the 
sphere was controlled by the circulation channel outside the phantom box 1210, which 
included the pump 1250 and the be^er 1240 for ink titration. In the experiment, a slab 
geometry was formed, with the slab having a thickness of 5.1 1 cm. The center of the ball 
was located 2.35 cm from the bottom of the slab. There were sources on the top plate (z 
= 0 cm) and detectors on the bottom plate (z = 5.1 1 cm). The signal to the sources was 
delivered via the source fibers 1270. The signal was light with a wavelength of 830 nm 
and was modulated at 70 MHz. The amphtude and phase of the diffusely transmitted 
hght were measured. Data was collected firom the detectors via the detectOr fibers 1220 
into the radio frequency imager 1230. The processing arrangement 1290 processed the 
results and generated the images. 

[0107] The images were fust generated using the Rytov approximation both 

excluding and including the reconstruction of the soiu-ce and detector coupling factors. 
Images were median filtered using standard techniques because this was found effective 
in removing speckle artifacts. FIG. 13A shows the images obtained excluding the 
reconstruction of the source and detector coupling factors, and FIG. 13B shows the 
images generated including the reconstruction. The X and Y ranges are in centimeters, 
and Z-shces are vertically arranged from 0.01 (top left) to 5.501 cm (bottom right) in 
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0.50-cin steps. It is observed that including reconstruction of the source and detector 
coupling factors improves image quality. In particular, the detail of the heterogeneity 
sphere can be seen in FIG. 13B for Z=2.01 to Z=3.51 cm, but it may not be observable 
over the corresponding Z-range in FIG. 13 A, 

[0108] It is to be understood that while the invention has been described in 

conjunction with the detailed description thereof, the foregoing description is intended to 
illush-ate and not limit the SCOpC Of the invention, which is defmed by the scope of the 
appended claims. Other aspects, advantages, and modifications are within the scope of 
the following claims. 
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